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ABSTRACT 

The  successful  application  of  holographic  interferometry  to  the  study 
of  three-dimensional  density  fields  around  bodies  in  wind  tunnel  experi- 
ments has  been  reported  in  the  literature  along  with  the  associated  math- 
ematical reduction  processes  for  the  basic  interf erometric  equation. 

The  present  report  has  extended  the  application  of  holography  as  a 
flow  visualization  technique  by  investigating  Mach  2.8  flow  over  a  rear- 
facing  l/8  in.  step.  Two  single- exposed  holograms  of  the  model  and  two 
double- exposed  frozen  fringe  holographic  interf erograms  of  the  flow  were 
created. 

Fast  Fourier  Transform  (FFT)  methods  have  been  reported  as  signifi- 
cantly reducing  computational  time  to  obtain  discrete  and  inverse  dis- 
crete Fourier  transforms.  Two  computer  programs  were  designed  to  apply 
FFT  methods  to  the  Fourier  transform  approach  to  the  basic  interf eromet- 
ric equation.  These  programs  were  used  to  investigate  an  alternate  way 
to  calculate  the  value  of  the  inside  integral.  Total  agreement  was  not 
found,  and  further  analysis  as  to  the  accuracy  of  each  method  is  needed. 
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I.   INTRODUCTION 

This  thesis  contains  "both  an  experimental  section  and  a  theoretical 
section.  The  experimental  portion  involved  the  construction  of  a  modified 
Z-type  Schlieren  system  to  conduct  holographic  investigation  of  a  rear- 
facing  step  in  supersonic  flow.  The  theoretical  portion  involved  inves- 
tigating current  data  reduction  techniques  for  determining  the  three- 
dimensional  distribution  of  a  field  "by  applying  fast  Fourier  transform 
methods  to  solve  the  basic  interferometric  equation.  Two  original  com- 
puter programs  were  devised  in  an  attempt  to  reduce  computer  computation- 
al time  required  for  interferometric  data  reduction. 


II.   INTRODUCTION  TO  HOLOGRAPHY 

Holography  is  a  three-dimensional  photography-like  process  which 
involves  the  recording  of  both  phase  and  amplitude  information  of  light 
waves  diffracted  by  an  object  illuminated  with  coherent  light.   Original- 
ly proposed  by  Gabor  £l,2]  in  19kQ,   holography  awaited  the  development  of 
the  laser  around  i960  to  become  an  explosive  field  of  its  own.  The  laser 
contributed  a  basically  monochromatic,  coherent  light  source.  The  prop- 
erty of  coherence,  or  the  capability  of  interference  with  itself,  enables 
the  recording  of  both  phase  and  amplitude  information. 

Photography  records  the  intensity  of  illumination,  which  is  the  time 
average  of  the  square  of  the  amplitude  of  the  light  field  at  any  point. 
In  holography,  a  coherent  light  beam  is  separated  into  a  reference  beam 
and  an  object  (or  scene)  beam.  The  reference  beam  is  directed  to  the 
holographic  plane  via  an  optics  system  while  the  object  beam  reaches  the 
same  point  after  being  diffracted  by  a  subject  of  interest.  The  emulsion 
on  the  holographic  plate  records  the  intensities  of  the  illuminating 
beams  and  the  interference  between  the  two  beams.  The  interference  com- 
ponent contains  the  phase  information  which  permits  reconstruction  of  an 
an  image  possessing  full  three-dimensionality  and  parallax  properties. 
No  attempt  will  be  made  to  present  the  full  mathematical  and  historical 
development  of  holography  which  is  amply  done  in  many  excellent  texts.  [3] 

Among  its  many  engineering  applications,  holography  is  well  suited 
for  flow  visualization.   Recent  successful  applications  of  holography  at 
the  Naval  Postgraduate  School  include  the  study  of  a  free  jet  and  the  de- 
termination of  the  density  field  about  a  cone  tilted  in  supersonic  flow 
in  a  conventional  wind  tunnel.  The  use  of  double-pulsed  holography  is 
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especially  helpful  as  a  tool  for  measurement  in  high  velocity  flows  near 
vibrating  equipment  such  as  free  jets  and  blow-down  wind  tunnels.  A  Q- 
switch  ruby  laser  is  used  as  a  high  power  energy  source  emitting  a  coher- 
ent light  pulse  on  the  order  of  several  joules  in  an  extremely  small  time 
span,  on  the  order  of  twenty  nanoseconds.  Vibrations  of  the  object  or  of 
the  optical  system  components,  which  must  remain  stable  within  wave  lengths 
of  light  during  exposure  of  the  hologram,  are  of  negligible  importance  due 
to  the  very  short  duration  of  the  exposure  time. 

A  double-pulsed  (or  double- exposed)  approach  is  advantageous  in  flow 
visualization.  In  the  Mach-Zender  interferometer,  comparison  of  the  ref- 
erence and  object  beams  are  made  at  the  same  time.   With  holography  it  is 
possible  to  separate  the  simultaneous  comparison  in  time.  Two  holograms, 
or  exposures,  are  made  on  the  same  holographic  plate  without  intermediate 
development.  The  first  exposure  is  made  with  no  flow  and  the  second  with 
flow.  The  individual  holograms  interfere  with  each  other,  creating  an 
interf erogram  of  the  flow.   Essentially,  this  interferogram  contains  only 
the  differences  in  optical  path  lengths  between  the  two  exposures.  There- 
fore, any  imperfections  in  optical  components,  (which  are  identical  for 
both  exposures),  are  subtracted  out,  thereby  eliminating  the  need  for  high 
quality  optics  \jf\. 

Another  method  of  taking  finite  fringe  holograms  is  the  double- exposure 
option  of  the  dual  hologram  method,  which  involves  using  two  holographic 
plates.   One  plate  is  used  for  each  exposure.   A  dual  hologram  holder  po- 
sitions each  plate  in  the  same  location  for  its  exposure.  The  improvement 
comes  from  the  versatility  of  almost  an  infinite  set  of  reconstruction 
possibilities  from  one  set  of  plates.   A  dual  hologram  holder  has  recently 
been  constructed  at  the  Naval  Postgraduate  School  to  demonstrate  its  ap- 
plication.  This  topic  was  the  subject  of  AIAA  Paper  73-210  [5]. 
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III.   INVESTIGATION  OF  A  REAR-FACING  STEP 

Holographic  investigation  of  a  rear-facing  step  in  supersonic  flow 
was  selected  for  two  main  reasons.  First,  extensive  theoretical  and  ex- 
perimental investigations  of  a  rear-facing  step  in  two-dimensional  super- 
sonic flow  are  available.  This  provides  ample  information  for  verifying 
results  from  this  experiment.  Secondly,  previous  experimental  results  do 
not  completely  agree  with  theoretical  descriptions.   A  source  of  conflict 
is  the  lip  shock*  s  influence  on  the  rest  of  the  flow  as  it  pertains  to  the 
reattachment  point  and  to  the  tail  shock.   With  this  in  mind,  it  was  hoped 
to  provide  a  hologram  which  could  be  used  to  accurately  reproduce  the  flow 
field. 

The  classical  flow  model  for  supersonic  base  flow,  (rear-facing  step 
flow)  ,  was  proposed  by  Chapman  \_6~]   in  1957  and  independently  by  Korst  [7] 
for  turbulent  flow,  and  is  commonly  known  as  the  Chapman-Korst  flow  model. 
Refer  to  Figure  1  for  a  sketch  of  the  Chapman-Korst  flow  model.   Roache 
[_ Q]   presented  a  literature  search  and  a  review  of  discrepancies  with  this 
model  in  his  doctoral  thesis  for  the  University  of  Michigan  in  I968.  The 
main,  unresolved  discrepancies  in  the  Chapman-Korst  flow  model  were:   that 
the  experimental  pressure  at  recompression  was  not  equal  to  the  theoreti- 
cal pressure  of  Region  Four  (see  Fig.  l);  and  that  the  isentropic  turn 
from  Regions  One  to  Two  (see  Fig.  l)  and  the  constant  pressure  assumptions 
were  invalidated  by  the  presence  of  the  lip  shock  (see  Fig.  2)  at  higher 
Mach  numbers.  L9j 

The  opportunity  to  provide  some  accurate  data  on  these  unresolved 
areas  of  interest  appeared  to  be  a  worthwhile  application  of  holography. 
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IV.   EXPERIMENTAL  MODEL 

The  actual  rear-facing  step  model  used  was  limited  in  size,  with  the 
maximum  step  size  being  one-eighth  inch.  Due  to  the  small  step  size,  the 
boundary  layer  thickness  (fi)  was  compared  with  the  step  size  (t)  to  de- 
termine if  distortions  from  the  boundary  layer  would  invalidate  experi- 
mental results.  From  Schlicting  [lo],  the  applicable  equations  are: 
6  =  5  (v  L  /  U)T  for  laminar  flow;  and  6  /  L  =  .3?  (  U  L/  v)"1'-5  for  the 
boundary  layer  thickness  at  the  trailing  edge  of  a  flat  plate  at  zero  in- 
cidence in  parallel  turbulent  flow.  The  6  /  t  ratio  results  were  .0538 
for  laminar  flow  and  .218  for  turbulent  flow.  It  was  felt  that  the  ac- 
tual flow  situation  was  laminar  and  that  this  step  size  would  still  offer 
meaningful  results. 

A  photograph  of  the  model  is  included  as  Figure  3»  Refer  to  Figure  h 
for  a  sketch  showing  the  dimensions  of  the  model. 
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V.   DESCRIPTION  OF  APPARATUS 

Figure  5  shows  a  complete  sketch  of  the  experimental  setup.   A  modified 
Z-type  Schlieren  system,  it  was  constructed  specifically  for  this  thesis. 
A  description  of  the  integrated  setup  will  be  followed  "by  listing  of  im- 
portant supportive  equipment. 

A  high  energy  pulse  was  discharged  from  a  Q-switch  ruby  laser  and 
passed  through  a  diffusing  lens  to  a  collimating  mirror  in  order  to  expand 
the  light  beam  to  sufficiently  illuminate  the  test  section.  The  colli- 
mated  beam  reflected  off  a  flat  mirror  towards  the  test  section.  Just 
prior  to  the  test  section,  the  beam  was  divided  by  a  beam  splitter  mount- 
ed on  a  tripod.  The  object  (scene)  beam  transmitted  through  the  beam 
splitter,  through  the  test  section,  and  reflected  off  another  flat  mirror 
towards  a  converging  mirror.  The  object  beam  was  then  reflected  towards 
the  holographic  plane  which  was  positioned  in  the  diverging  part  of  the 
beam  where  the  test  object  was  in  focus.  The  part  of  the  collimated  beam 
reflected  by  the  beam  splitter,  the  reference  beam,  reflected  off  three 
flat  mirrors  and  also  illuminated  the  holographic  plane. 

Other  equipment  included  an  air  supply  for  the  tunnel  provided  by  an 
air  compressor;  a  power  supply  and  cooling  system  for  the  laser;  a  15  mW 
continuous  wave  Helium-Neon  laser  used  for  alignment  of  the  Q-switched 
ruby  laser  with  the  optics  system;  an  autocollimator  for  internal  align- 
ment of  the  Q-switch  laser;  a  movable  platform  for  the  laser  so  it  can  be 
elevated  to  the  proper  height  for  the  optics  system;  a  table  for  the  holo- 
gram holder;  and  mountings  and  supports  for  all  optic  components.  Figures 
6-9  are  pictures  showing  the  experimental  setup  from  four  different  views. 

Safety  equipment  included  eye  protectors  and  sound  attenuators. 
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VI.      EXPERIMENTAL  PROCEDURE 

Appendix  A  contains  a  detailed  description  of  the  entire  experimental 
procedure  from  alignment  of  system  components  to  exposure  of  various  types 
of  holograms. 

There  are  two  holographic  methods  that  can  be  applied  to  flow  visual- 
ization in  wind  tunnels.  The  first  method  is  commonly  referred  to  as  the 
real  fringe  method,  and  the  second  as  the  frozen  fringe  or  double- exposed 
method. 

Real  fringe  methods  require  exposure  of  the  holographic  plate  with  no 
flow  in  the  wind  tunnel  and  either  development  of  the  hologram  in  place 
or  precise  repositioning  of  the  developed  plate.      As  flow  is  increased  in 
the  wind  tunnel,   an  evolution  of  real  fringe  patterns  could  be  observed 
through  the  hologram.      This  method  was  not  used  in  this  experiment  due  to 
the  unavailability  of  in-place  developing  equipment  or  a  dual-hologram 
holder. 

Double- exposed  or  frozen  fringe  methods  require  two  exposures  of  the 
same  holographic  plate.-     The  first  exposure  is  taken  with  no  flow  in  the 
wind  tunnel,    and  the  second  with  flow.     The  interference  between  the  two 
exposures  results  in  a  double-pulsed  (exposed)   interf erogram.      In  this 
experiment  both  exposures  were  made  on  the  same  plate  with  the  second  ex- 
posure occurring  while  Mach  2.8  flow  existed  in  the  test  section. 
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VII.      FILM  AND  DEVELOPMENT  PROCEDURE 

There  are  essentially  two  sources  of  holographic  films:     Eastman  Kodak 
Company  and  Agfa-Gevaert.      Agfa-Gevaert  films  were  used  exclusively  in 
this  experiment  because  their  75  series  film  is  designed  for  maximum  sen- 
sitivity around  the  emission  spectrum  of  ruby  lasers  (69^3  A).      Due  to 
the  short  duration  of  the  exposure  time  of  the  Q-switch  ruby  laser,   about 
10-20  nanoseconds,    10E75  film  was  used  because  8E75  film  turned  out  to  be 
too  slow  in  recording  information. 

Development  of  holograms  can  include  various  features  such  as  bleach- 
ing, prehardening,  or  varying  development  time.     The  various  alternatives 
are  examined  in  Ref.    11.     For  this  experiment  the  following  development 
procedure  was  used: 

1.  Five  minutes  in  Kodak  D-19  Developer 

2.  Thirty  seconds  in  an  acetic  acid  stop  bath  of  standard  dilution. 

3.  Five  minutes  in  a  standard  fixer. 

4.  Five  minutes  of  washing  followed  by  immersion  in  a  wetting  agent 
(Kodak  Photoflo)  prior  to  drying. 
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VIII.   RECONSTRUCTION  TECHNIQUE,  EQUIPMENT,  AND  RESULTS 

A  15  mW  continuous  wave  Helium-Neon  laser,  previously  used  for 
alignment,  was  the  coherent  light  source  used  for  reconstruction.  A  col- 
limating  lens  was  attached  to  this  laser  in  order  to  expand  the  recon- 
struction "beam  to  approximately  two  inches.   Other  reconstruction  equip- 
ment included  a  hologram  holder,  type  52P  Polaroid  film,  a  graphic  view  ' 
camera,  a  tripod,  and  a  black  drape. 

The  two-inch  diameter  collimated  beam  incident  to  the  hologram  created 
two  images,  whether  projected  images  are  real  or  virtual  is  discussed  in 
Chapter  3  of  Ref.  3  for  the  various  combinations  of  converging,  diverging, 
and  collimated  reference  and  scene  beams.  In  this  experiment,  approximate 
symmetry  about  the  vertical  axis  of  the  holographic  plate  was  attempted. 
This  resulted  in  the  real  image  being  projected  as  a  converging  wave  front 
and  in  the  virtual  image  being  projected  as  a  diverging  wave  front. 

In  order  to  focus  the  real  image,  the  lens  was  removed  from  a  camera 
and  the  image  was  projected  on  the  ground  glass  backing  of  the  camera. 
Pictures  of  these  projections  were  taken  at  exposure  times  which  varied 
with  the  intensity  of  the  real  image  beam.   Normal  exposure  times  were 
approximately  one  second. 

Figure  10  is  a  sketch  of  the  reconstruction  setup  and  Figures  11-14 
are  examples  of  photographs  of  the  projections  of  the  real  image  of  holo- 
grams taken  with  this  experimental  setup. 

Figures  13  and  14  represent  the  final  experimental  attempts  to  create 
a  frozen  fringe  hologram  of  sufficient  quality  to  warrant  data  reduction. 
The  edges  of  the  experimental  model  show  that  both  holograms  are  slightly 
out  of  focus.  These  holograms  were,  however,  sufficient  to  demonstrate 
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the  applicability  of  holography  to  flow  visualization  of  a  rear-facing 
step  in  supersonic  flow.     Additionally,   the  shadowgraph  resulting  on  the 
developed  plate  showed  the  basic  features  of  the  Chapman-Korst  flow  model. 
Sufficient  information  was  not  obtained  to  analyze  any  of  previously  men- 
tioned discrepancies  of  that  model. 

The  remaining  effort  of  this  thesis  focused  on  the  possibilities  of 
improving  current  data  reduction  methods. 
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IX.   CURRENT  DATA  REDUCTION  TECHNIQUES 

There  are  currently  several  approaches  for  reducing  an  asymmetric 
density  field.  One  method  involves  a  series  expansion  of  a  complete  set 
of  orthogonal  polynomials  which  are  invariant  to  a  rotation  of  coordinates 
£12^.  Matulka,  [.13]  used  this  method  and  created  a  computer  program  which 
took  approximately  two  hours  of  computer  time  to  reduce  the  data  for  the 
asymmetric  case. 

A  second  approach  uses  Fourier  transform  techniques  suggested  by 
Rowley  [_lk]   and  further  developed  by  Junginger  and  van  Haeringen  |~15]« 
This  approach  is  thoroughly  discussed  by  Sweeney  [l6]  in  his  doctoral  the- 
sis for  the  University  of  Michigan  in  1972. 

Closely  following  Sweeney's  development  and  nomenclature,  the  funda- 
mental integral  equation  to  be  solved  in  three-dimensional  interferometry 
is 

/  f(x,y,z)  dS  =  f  (9-1) 

Si 

where  f(x,y,z)=n(x,y,z)-n  ;  f(x,y,z)  is  the  refractive  index  relative  to 
a  known  value  n  ,  S.  is  the  distance  along  the  i   ray  traversing  the  test 
section,  f>.    is  the  known  path  length  function,  and  the  integral  is  eval- 
uated between  boundaries  imposed  by  experimental  configuration  or  by  prior 
knowledge  regarding  f(x,y,z).  The  objective  is  to  reconstruct  a  real 
scalar  function,  the  relative  refractive  index  field,  from  the  scalar  val- 
ues of  line  integrals  through  the  field.   When  the  gradient  of  the  refrac- 
tive index  variation  is  small  enough  to  produce  negligible  ray  curvature, 
eqn.  (9-l)  alone  is  the  basic  interf  eromet.^1  c   equation  to  be  solved.  The 
line  integrals  are  then  evaluated  along  straight,  not  curved,  lines. 
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The  presentation  of  the  problem  is  simplified  by  the  suppression  of 
one  of  the  dimensions.     Solutions  for  the  two-dimensional  case  can  be  lay- 
ered to  reproduce  the  three-dimensional  distribution,   thereby  permitting 
simplification  without  loss  of  information. 

The  use  of  polar  coordinates  on  the  two-dimensional  equivalent  of 
eqn.    (9-1)   yields 

J"  f(x,y)   dS  =  p(/?,e)  (9-2) 

where  p  varies  to  cover  the  domain  of  f(x,y)   and  0  varies  between  ±  rr/2. 
Eqn.    (9-2)    then  becomes  the  fundamental  integral   equation  to  be  solved. 

If  the  refractive  index  field  is  radially  symmetric  or  if  it  is  inde- 
pendent of  one  of  the  coordinates,   then  the  solution  of  eqn.    (9-2)   is 
known  and  has  been  extensively  used  in  interferometric  studies.      For 
f(x,y)=f(y)    only,    eqn.    (9-2)   becomes 

L/2 

V  (y)   =  S         n(y)    dx  =  L  n(y)  (9-3) 

'  -L/2 

and  the  solution  is  n(y)   =  ®(y)  L. 

For  the  radially  symmetric  case,   only  one  field  of  view  is  necessary 

and  eqn.    (9-2)   becomes 

*(x)   =  2     r f(r)   dy  =  2  T  ^J*  (9-4) 

/  r=x  x     (r  -x  )2 

which  is  a  form  of  Abel's  integral  equation. 

The  solution  of  eqn.    (9-4)  may  be  obtained  by  converting  the  integral 

into  a  convolution  integral  by  using  appropriate  transformations.     The 

solution,   f (r) ,    then  becomes 

*»-$  z'V-V    s[#-]"  CM) 

Sweeney  uses  the  definitions  of  Fourier  transforms,  a  few  fundamental 
transform  pairs,  the  procedure  for  inverting  the  Fourier  transformed 
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equivalent  of  eqn.    (9-2)    as  outlined  by  Junginger  and  van  Haeringen  [l5]> 
and  an  introduction  of  polar  coordinates  to  the  transformed  equation  to 
yield 


where  "*"  denotes  the  convolution.     Expanding  the  convolution  integral 
and  replacing  z  with  r  sin(a>-0),   the  final  result  for  f(r,»)   is: 

f(r,p)  .  ^      f  d9    j*-  Ufefe^V  (9-7) 

2tt     -n/2  -oo  v      '     r 

This  is  the  basic  relationship  between  the  pathlength  function  and  the 
refractive  index  field. 

Referring  to  eqn.  (9-6) ,  the  integral  over  0  is  a  linear  summation  of 
the  contributions  for  each  direction  of  view.  The  contribution  from  each 
value  of  0  is  equal  to  the  convolution  of  the  derivative  of  the  path 
length  function  and,  (l/z)  with  (l/z)  evaluated  at  z  =  r  sin(cp-Q) 

Wheny=  r  sin(<p-0)  ,  the  integrand  of  eqn.  (9-7)  becomes  singular. 
This  presents  no  problem  since  the  contributions  to  the  integral  near  the 
origin  tend  to  cancel  for  an  odd  function  (like  l/z).  Thus  the  integral 
remains  finite. 

The  integration  over  p  in  eqn.  (9-7)  is  a  Hilbert  transform  of  ®(o,q) 
and  the  properties  of  Hilbert  transforms  along  with  the  use  of  discrete 
data,  guarantee  the  boundness  of  eqn.  (9-7).  Both  this  form  of  the  inte- 
gration over  p  and  Van  Houten's  [_17~]   modification  to  it  will  be  called 
the  inside  integral  in  the  next  section. 

A  finite  number  of  values  of  the  frequency  plane  must  be  determined. 
One  method  for  their  determination  uses  linear  interpolation  from  adja- 
cent radial  lines  of  a  superimposed  grid  system.   Success  of  this  method 
depends  on  the  accuracy  of  the  interpolation.   If  the  pathlength  data  is 


21 


taken  for  a  sufficiently  large  number  of  angular  orientations,  linear 
interpolation  may  be  used. 

These  linearly  interpolated  frequency  plane  values  (f{ptQ))   must  be 
convoluted  with  l/z,  i.e. 

[?(/.e)  *  l/z]  (9-8) 

with  "*"  meaning  convolution. 

Using  the  following  notation  of  F  for  obtaining  the  discrete  Fourier 

•i 
transform  (IFT)  and  F   for  obtaining  the  inverse  DFT  or  I EFT,  the  re- 
maining steps  required  to  perform  the  integration  over p   are  summarized 
as  follows: 

The  DFT  of  the  convolution  is  obtained 

Fty(p,e)  *  l/z]  (9-9) 

followed  by  the  differentiation  of  this  function 

d(F[y(/.6)  *  l/z])  (9-10) 

Finally  the  value  for  the  integration  over  p  is  obtained  by  taking  the 
IIFT  of  the  preceeding  function. 

F"1  [d(F[p(^,e)  *  l/z])]  (9-11) 

Van  Houten  \Vf\   used  the  linear  interpolation  approach  to  obtain  the 
frequency  plane  values  by  a  computer  program  designed  to  calculate  the 
three-dimensional  density  field  from  data  obtained  from  holographic  in- 
terferometry.  The  frequency  plane  values  for  the  density  field  are  rep- 
resented by  expanded  fringe  shift  information  calculated  from  the  fringe 
number  of  specific  points  on  the  interferogram.  Van  Houten  modified  eqn. 
(9-7) i  the  basic  interf erometric  equation  to  the  following: 

_±  3RMAX     z(po+pr)  +  dpo-p,<)-2r,(Fo,r) 

f(x0»yJ  -  — ?   *  L  s        5 do 

2tt        0         0  c  ' 

2g(/>_  ,p) 

+  -3HMA3H   df  (9-12) 
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with  g(p,p>)   being  the  expanded  fringe  shift  information  and  the  integral 
inside  the  "brackets  referred  to  as  the  inside  integral  instead  of  the  in- 
tegration over  p. 

Van  Houten' s  computer  program  uses  Cote's  sixth  order  method  Ll8^  for 
integrating  the  first  term  of  the  inside  integral.  He  applied  a  correct- 
ing term  for  the  contribution  of  the  integral  from  3RMAX  to  infinity. 
The  sum  of  these  two  terms  is  called  the  inside  integral.  If  NP  is  the 
number  of  data  points  for  each  view,  and  if  K  is  the  number  of  views, 

the  inside  integral  must  be  evaluated  (NP  •  K)  times  with  Van  Houten' s 

2 

computer  program.   Additionally,  for  each  of  the  NP  grid  intersections, 

K  interpolations  for  G(y?  ,©,  )  are  required.  The  subsection  of  Van  Houten' s 
computer  program  which  calculates  the  value  of  the  inside  integral  is  in- 
cluded in  Appendix  B.  The  next  section  will  investigate  alternatives  for 
this  part  of  Van  Houten' s  computer  program. 
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X.   APPLICATION  OF  THE  FAST  FOURIER  TRANSFORM 

Fast  Fourier  transform  (FFT)  methods  are  simply  a  more  efficient  means 
of  obtaining  the  discrete  Fourier  transform  (DFT).  Their  increased  effi- 
ciency is  due  to  the  significant  reduction  in  the  number  of  complex  math- 
ematical operations  required.  For  a  function  of  N  points,  normal  calcula- 

2 

tion  of  the  DFT  requires  nearly  N     operations  while  the  FFT  radix-2  method 

reduces  this  to  (N/2)   log?(N)   complex  multiplications,    (n/2)   log„(N)   com- 
plex additions,   and  (N/2)   log?(N)    complex  subtractions.      As  can  be  seen 
from  Figure  15,   for  N=1024,   this  results  in  a  computational  savings  of 
more  than  200  to  1.   [19]     Higher  radix  FFT  methods  further  reduce  the  num- 
ber of  operations  required. 

Operations  usually  associated  with  FFT  applications  include:      (l)    com- 
puting a  spectrogram;    (2)    computing  the  convolution  of  two  time  series; 
and  (3)   computing  the  correlation  of  two  time  series.     The  second  of  these 
operations  has  been  shown  to  have  direct  application  to  the  calculation 
of  the  value  of  the  inside  integral  for  determining  the  three-dimensional 
density  field.     Reference  20  gives  a  thorough  discussion  of  the  FFT  and 
the  required  steps  for  convolution  and  correlation.     These  detailed  pro- 
cedures for  computing  the  convolution  are  included  in  Appendix  C  for  the 
following  combination  of  functions:      (l)   two  finite  length  functions;    (2) 
a  finite  and  an  infinite  length  function  using  the  select-savings  method; 
and  (3)    a  finite  and  an  infinite  length  function  using  the  overlap-add 
method. 

Two  original  Fortran  programs  were  devised  to  apply  FFT  methods  to  the 
calculation  of  the  values  of  the  inside  integral.      Specifically,   alterna- 
tives to  Mode  One  operation  of  Van  Houten's  computer  program  for  the 

zh  1 


2  2 

asymmetric  Gaussian  case  of  f(z)  =  (l-y  )    exp  (-2x  )   for  65  points  from 

-1.5  <  R  <  1.5  were  investigated.      Mode  One  operation  of  Van  Houten's 
program  requires  refractive  index  information  supplied  "by  SUBROUTINE  IN- 
PUT.    Using  a  finite  sample  of  65  points  with  an  interval  between  samples 
of  . 0W>&75 »  Van  Houten's  program  required  1.3388  seconds  of  computer  time 
just  to  calculate  the  value  of  the  inside  integral.     This  time  was  ob- 
tained by  using  computer  library  subroutines  SETIME  and  GETIME  placed 
just  prior  to  and  immediately  after  the  subsection  of  Van  Houten' s  com- 
puter program  which  calculates  the  value  of  the  inside  integral. 

Van  Houten's  program  also  uses  a  SUBROUTINE  COEFF  to  calculate  coeffi- 
cients used  in  various  parts  of  the  program.     M0DE1 ,   a  modification  of  Van 
Houten's  program  which  eliminated  the  need  for  SUBROUTINE  COEFF  by  return- 
ing its  calculations  to  the  main  program,   is  included  in  Appendix  D.      This 
modification  includes  the  essentials  of  Van  Houten's  thesis  with  the  fol- 
lowing corrections:      Gl  =   .292857  instead  of   .2292857;    several  divisions 
by  59.296  were  corrected  to  57.296;    and  all  plotting  directions  were  de- 
leted due  to  the  storage  space  limitation  of  the  IBM  36O  on-line  time 
sharing  capability.     This  modification  followed  the  work  done  by  Everett 

[21]. 

Both  FFT  method  programs  used  a  finite  sample  of  the  expanded  fringe 
shift  information  obtained  from  Van  Houten's  computer  program.  His  meth- 
od for  obtaining  these  frequency  plane  values  followed  the  linear  inter- 
polation approach  as  described  by  Sweeney.   Sixty-five  finite  sample 
points  were  obtained  and  they  had  to  be  convoluted  with  the  infinite 
function  l/z.   J0NES1  numerically  convoluted  these  functions  while  J0NES2 
used  the  FFT  select-savings  method  to  obtain  the  convolution.   Once  the 
convolution  was  obtained,  both  programs  followed  the  procedures  described 
in  the  preceding  section,  except  that  FFT  methods  were  used  to  obtain  all 
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DFT's  and  IDFT's.  The  two  separate  approaches  for  calculating  the 
convolution  served  as  an  internal  check  that  the  correct  convolution  was 
being  obtained.  The  complete  programs  are  included  in  Appendix  F.  along 
with  a  listing  of  the  expanded  fringe  shift  information  obtained  from  Van 
Houten's  program  and  used  as  frequency  plane  values  for  both  programs. 

The  numerical  methods  used  in  J0NES1  required  multiplying  the  final 
result  by  the  appropriate  weight  function,  the  interval  between  sampled 
points  of  the  expanded  fringe  shift  information.  J0NES2  used  the  FFT 
select- savings  method  to  calculate  the  convolution.  The  dimension  of  the 
l/R  vector  was  chosen  to  be  512  following  Brigham's  [20^  recommendation 
as  the  optimal  choice  for  values  of  Q  between  50  and  99.  His  table  of 
recommendations  for  the  size  of  arrays  needed  for  FFT  convolution  is  in- 
cluded in  Figure  16.  The  sample  interval  for  both  of  these  FFT  programs 
was  selected  to  agree  with  the  interval  used  in  Van  Houten's  program, 
namely  .046875.  For  N  =  512,  and  DR  being  defined  as  2RMAX/N,  an  R  of 
12.0  was  required.  Using  FFT  methods  to  obtain  the  Fourier  coefficients 
of  the  convolution  also  required  a  weighting  function.  In  the  complex 
multiplication  of  the  normalized  Fourier  coefficients  of  the  expanded 
fringe  shift  information  and  the  l/R  vector,  the  resulting  Fourier  coef- 
ficients incorporated  an  extra  division  by  N.  Instead  of  a  weighting 
function  of  DR  then,  2R  was  chosen  as  the  appropriate  weighting  function. 
When  the  IDFT  of  the  Fourier  coefficients  was  obtained,  the  values  of  the 
convolution  by  this  FFT  method  were  in  exact  agreement  with  the  convolu- 
tion obtained  by  numerical  methods. 

The  differentiation  technique  used  in  J0NES1  and  J0NES2  required  an 
increment  between  steps  of  differentiation  that  was  equal  to  the  sample 
data  interval.   Thus,  the  same  DR  was  used  for  this  process. 
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All  FFT  calculations  in  both  programs  utilized  SUBROUTINE  CFFT2,   an 
algorithm  for  computing  the  mixed  radix  fast  Fourier  transform  [22^.     The 
complete  single  precision  version  of  CFFT2  is  included  in  Appendix  E. 
Chosen  for  its  versatility,   it  contains  a  detailed  and  straight-forward 
usage  description  at  the  "beginning  of  the  subroutine.     Several  known  Fou- 
rier transform  pairs  were  investigated  to  ensure  that  CFFT2  was  being  used 
correctly.      Correct  DFT's  and  ITFT's  were  easily  obtained. 

As  was  previously  mentioned,  FFT  methods  involve  complex  mathematical 
operations.      When  using  only  real  functions,   as  was  the  case  with  the  ex- 
panded fringe  shift  samples  and  the  l/z  function,   the  frequency  plane 
values  could  be  placed  in  the  real  coefficient  array  of  CFFT2  and  the 
other  function  could  be  placed  in  the  imaginary  coefficient  array  of  the 
same  call  to  CFFT2.     This  would  result  in  the  significant  reduction  of 
time  and  storage  space  required  to  calculate  DFT's  and  IDFT's.     This  in- 
creased efficiency  was  not  utilized  in  either  J0NES1  or  in  J0NES2.     In- 
stead,  both  programs  are  quite  inefficient,    with  redundant  storage  used 
for  recording  and  checking  each  step  in  the  entire  process  of  calculating 
the  value  of  the  inside  integral.      If  one  only  desired  a  print-out  of  the 
inside  integral  values,    these  programs  could  easily  be  reduced  to  four 
matrices,    and  possibly  to  two  matrices  since  both  functions  were  real  in 
this  case.      Both  programs,   designed  as  possible  alternatives  for  Van 
Houten's  subsection  that  calculated  the  value  of  the  inside  integral, 
were  only  written  to  replace  inside  integral   calculations  for  Mode  One 
operation.     Hence,   in  addition  to  condensing  these  programs  into  effi- 
cient packages,   investigation  of  all  other  modes  of  operation  is  required 
before  a  complete  integration  of  FFT  methods  into  Van  Houten's  program 
can  be  achieved. 
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XL,      RESULTS  AND  RECOMMENDATIONS  FOR  FFT  APPLICATION 

The  computer  time  required  for  just  the  calculation  of  the  inside 
integral  subsection  of  Van  Houten's  program  for  Mode  One  operation  of  the 
asymmetric  Gaussian  case  took  1.3388  seconds.     Although  the  increased  ef- 
ficiency of  FFT  methods  are  significant  for  large  values  of  N,   it  was  not 
known  if  FFT  methods  would  reduce  computational  time  required  in  this 
case  where  N  was  only  64  sample  intervals.     Nevertheless,   if  time  had 
permitted  integration  of  a  condensed  and  efficient  form  of  either  J0NES1 
or  J0NES2  into  Van  Houten's  program,   it  was  felt  that  some  reduction  in 
computation  time  would  have  been  achieved. 

Both  approaches  to  obtain  the  convolution,   J0NES1  using  numerical 
methods  and  J0NES2  using  FFT  methods,   produced  convolution  values  in  ex- 
act agreement  with  each  other  and  inside  integral  values  within   .002  per 
cent  of  each  other.      A  listing  of  these  results  side  by  side  to  the  in- 
side integral  values  obtained  from  Van  Houten's  program  is  included  in 
Appendix  G.     There  was  general  agreement  between  the  values  of  the  inside 
integral  as  calculated  by  Van  Houten's  program  compared  with  those  ob- 
tained by  FFT  methods.      Near  a  value  of  R   equal  to  ±  .8  the  apparently 
significant  percentage  difference  between  the  two  methods  was  misleading 
since  the  magnitudes  involved  were  small.      However,   for  values  of  R  near 
the  discontinuity  in  the  density,    there  was  significant  disagreement  in 
the  values  of  the  inside  integral.      Further  analysis  as  to  the  accuracy 
of  the  methods  involved  is  required  to  resolve  this  problem. 

Shifting  the  expanded  fringe  shift  information  from  the  middle  of  its 
array  to  the  beginning  of  an  array  whose  remaining  points  were  filled 
with  zeros,    as  suggested  by  the  select-savings  method,   resulted  in  the 
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inside  integral  values  corresponding  to  the  -1.5  <  R  <  1»5  interval  being 
stored  in  the  output  array  beginning  with  the  (n/2  +  l)  position.  In 
this  case,  the  appropriate  values  of  the  inside  integral  were  stored  from 
the  257  through  the  321s  positions  of  the  output  vector.  According  to 
the  select- savings  method  (refer  to  Appendix  C.2),  the  first  "Q-2"  points 
of  the  output  vector  are  meaningless.  Q  is  defined  as  the  number  of  fre- 
quency plane  samples  used  in  the  calculation  of  the  convolution. 
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XII.      SUMMARY  OF  RESULTS 

An  experimental  setup  was  designed  and  constructed  to  conduct 
holographic  investigation  of  a  rear-facing  step  in  supersonic  flow.     Al- 
though the  double- exposed  (frozen)   holograms  produced  did  not  warrant 
further  analysis  of  discrepancies  of  the  Chapman-Korst  flow  model,   they 
were  of  sufficient  quality  to  demonstrate  an  application  of  holography 
to  rear-facing  step  flow  visualization  near  vihrating  equipment  such  as 
the  supersonic  "blowdown  wind  tunnel. 

Theoretical  investigations  of  current  data  reduction  techniques  for 
determining  the  three-dimensional  density  field  from  interf erometric  data 
led  to  investigation  of  the  applicability  of  fast  Fourier  transform  meth- 
ods to  reduce  computational  time  required  for  determining  the  inside  in- 
tegration over p  in  the  basic  interf erometric   equation.      Two  programs 
were  devised  which  calculated  the  value  of  the  inside  integral  using  FFT 
methods.      Although  the  convolution  was  obtained  differently  for  each  pro- 
gram,  numerically  in  J0NES1  and  with  FFT's  for  J0NES2,    the  convolution's 
were  identical.      Comparison  of  these  values  with  those  obtained  from  Van 
Houten's  program  lacked  close  agreement  for  the  ±  (1.07  -  1.31)   intervals 
for  R.      The  difference  between  these  methods  requires  further  analysis  as 
to  the  accuracy  of  the  methods.     Further  analysis  is  also  needed  to  effi- 
ciently reduce  JONES   1   and  J0NE32  and  to  integrate  the  programs  into  Van 
Houten's  program  for  all  of  its  modes  of  operation. 
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Figure  3.   The  Experimental  Model  of  a  Rear-Facing  Step 
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Figure  6.   A  view  of  the  left  side  of  the 
the  experimental  setup  showing  the  laser 
and  the  optics  system  up  to  the  beam 
splitter. 


Figure  7.   View  of  the  left  side  of  the  test 
section,  also  showing  the  Q-switched  laser, 
the  beam  splitter,  and  mirrors  which  direct 
the  reference  beam  over  the  test  section. 
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Figure  8. 

A  view  of  the  right  side 
of  the  test  section. 


Figure  9 


A  view  of  the  right  side 
of  the  experimental  setup 
showing  the  holographic 
plane  and  the  final  mirrors 
for  the  reference  and  the 
scene  beams. 
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Figure  11. 

Photograph  of  the  pro- 
jected real  image  of  a 
single-exposed  hologram. 


Figure  12. 

Photograph  of  the  pro- 
jected real  image  of  a 
single-exposed  hologram. 
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Figure  13. 

Photograph  of  the  pro- 
jected real  image  of  a 
double-exposed  hologram 
of  a  rear-facing  1/8  in. 
step  in  Mach  2.8  flow. 


Figure  14. 

Photograph  of  the  pro- 
jected real  image  of  a 
double-exposed  hologram 
of  a  rear-facing  1/8  in, 
step  in  Mach  2.8  flow. 
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Figure  15.  The  Number  of  Operations  Required 
to  Compute  the  Discrete  Fourier  Transform  by 
Direct  Calculation    Compared  with  Using   the 
FFT     Algorithm. 
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10  32  5 

11-19  64  6 

20-29  128  7 

30-49  256  8 

50-99  512  9 

100-199  1024  10 

200-299  2048  11 

300-599  4096  12 

600-999  8192  13 

1000-1999  16384  14 

2000-3999  32768  15 

Optimum  Value  of  N  for  FFT  Convolution, 

FIGURE  16. 
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APPENDIX  A.   DETAILED  EXPERIMENTAL  PROCEDURES 

I.  ALIGNMENT  OF  SYSTEM: 

a)  Turn  on  cooling  system  and  warm  up  power  supply. 

t>)  Check  internal  alignment  of  laser  with  autocollimator. 

c)  Insert  2.5  mm  orifice  inside  laser  cavity  in  order  to  eliminate 
higher  order  transverse  mode  emission. 

d)  After  ensuring  all  endangered  personnel  are  wearing  eye  protec- 
torsf  fire  the  Q-switch  laser  to  "burn  a  spot  on  underexposed  (black)  film. 

e)  Make  a  needle  hole  in  the  center  of  the  "burn  spot  and  fire  the  la- 
ser again  to  make  a  burn  spot  on  another  fixed  piece  of  underexposed  film. 
Always  check  for  personnel  wearing  eye  protectors  prior  to  firing  the  Q- 
switch  laser. 

f)  Align  a  15  mW  CW  laser  team  through  the  2.5  mm  orifice  and  through 
the  pin  holes  in  the  two  pieces  of  film  hy  using  a  flat  mirror. 

g)  Elevate  table  so  that  the  laser  beam  hits  the  first  mirror  in  the 
optics  system  in  its  center.  Ensure  that  the  focal  point  of  the  diffusing 
lens  is  located  at  the  focal  point  of  the  parabolic  mirror  in  order  to  ob- 
tain a  collimated  beam  entering  the  optics  system. 

h)  Align  optics  system,  ensuring  test  section  and  holographic  plane 
are  properly  illuminated. 

i)  Adjust  the  location  of  the  holographic  plane  to  ensure  that  the 
subject  image  is  in  focus  in  the  diverging  part  of  the  scene  beam. 

j)  Adjust  the  final  flat  mirror  in  the  reference  beam  in  order  to 
make  sure  path  lengths  of  the  two  beams  are  within  one-half  inch  of  each 
other. 
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k)  Remove  CW  laser  and  alignment  mirror  as  well  as  remaining  alignment 
pinholes. 

l)  Ready  to  fire. 

II.      EXPOSURE  OF  HOLOGRAPHIC  PLATES: 

a)  Single  Pulse  (Exposure): 

1)  Insert  plate  in  holder  underneath  the  black  drape. 

2)  Assistant  charges  the  power  supply  and  notifies  experimentor 
when  ready  to  fire. 

3)  Experimentor  lifts  the  drape  while  simultaneously  calling  to 
fire  the  laser,   and  replaces  the  drape  immediately  after  observing  lamp 
flash  of  the  Q-switch  laser. 

4)  Remove  the  exposed  holographic  plate  to  a  labeled  film  box 
while  still  under  the  drape  to  eliminate  further  exposure. 

5)  Hologram  is  ready  for  development. 

b)  Double-Pulse  (Exposure): 

1)  Repeat  steps  II. a. 1-3.  Wear  sound  attenuators. 

2)  Open  air  supply  and  create  desired  flow  in  the  test  section. 
(Mach  2.8  flow  in  this  experiment) 

3)  Repeat  steps  II. a. 2-5.  Assistant  shuts  off  air  supply  after 
firing  laser. 
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APPENDIX  C.      FFT  CONVOLUTION  PROCEDURES 

C.l     COMPUTATION  PROCEDURE  FOR  FFT  CONVOLUTION  OF  FINITE-LENGTH  FUNCTIONS. 

1".     Let  x(t)   and  h(t)   be  finite-length  functions  shifted  from  the 
origin  by  a  and  b,   respectively. 

2.  Shift  x(t)   and  h(t)   to  the  origin  and  sample: 

x(k)  =  x(kT  +  a)  k  =  0,1,... ,P  -1 

h(k)  =  h(kT  +  b)  k  =  0,1,. . . ,Q  -1 

3.  Choose  N  to  satisfy  the  relationships: 

N  >  P  +  Q  -  1 

N  =  2      integer  valued 
where  P  is  the  number  of  samples  defining  x(t)  and  Q  is  the  number  of  sam- 
ples defining  h(t)„ 

4.  Augment  with  zeros  the  sampled  functions  of  step  (2). 

x(k)  =0     k  =  P,P  +  1,...,  N  -  1 
h(k)  =0     k=Q,Q+l,...,N-l 

5.  Compute  the  discrete  transform  of  x(k)  and  h(k) . 

N-l 
X(n)  =  E  x(k)  exp(-j2iTnk/N) 
k=0 

N-l 
H(n)  =  E  h(k)  exp(-j2n  nk/N) 
k=0 

6.  Compute  the  product: 

Y(n)  =  X(n)H(n) 

7.  Compute  the  inverse  discrete  transform  using  the  forward  trans- 
form: 

N-l  , 
y(k)  =  E  ^  Y"(n)  exp(-j2n  nk/N) 
n=0 

k6  ( 


C.2     COMPUTATION  PROCEDURE  FOR  FFT  CONVOLUTION:      SELECT-SAVINGS  METHOD: 
i.     Let  Q  he  the  number  of  samples  representing  h(t). 

2.  Choose  N  according  to  the  table  in  Figure  16. 

3.  Form  the  sampled  periodic  function  h(k). 

h(k)  =  h(kT)  k  =  0,1 Q-  1 

=  0  k  =  Q,Q  +  1 N  -  1 

4.  Compute  the  discrete  Fourier  transform  of  h(k) 

N-l 
H(n)  =  E  h(k)  exp(-j2nnk/N) 
k=0 

5.  Form  the  sampled  periodic  function: 

x±(k)  =  x(kT)  k  =  0,1 N-l 

6.  Compute  the  discrete  Fourier  transform  of  x. (k) . 

N-l 
X.(n)  =  E  x.  (k)  exp(-j2Trnk/N) 
1     k=0 

7.  Compute  the  product: 

Yi(n)  =  Xi(n)H(n) 

8.  Compute  the  inverse  discrete  transform  of  Y. (n) . 

N-l 
yi(k)  =  E  7j  Y1(n)#exp(-j2nnk/N) 
n=0 

9.  Delete  samples  y. (0) ,y. (l) , . . .y. (Q  -  2),  and  save  the  remaining 
samples. 

10.  Repeat  steps  5  to  9  until  all  sections  are  computed. 

11.  Comhine  the  sectioned  results  by  the  relationships: 
y(k)  undefined         k  =  0,l,...,Q-2 

y(k)  =  yi(k)  k  =  Q  -  1.Q.....N  -  1 

y(k+N)  =  y2(k+Q-l)       k  =  0,l,...,N-Q 
y(k+2N)  =  y„(k+Q-l)       k  =  0,1,..., N  -  Q 

0 
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C.3     COMPUTATION  PROCEDURE  FOR  FFT  CONVOLUTION:      OVERLAP-ADD  METHOD. 

1.  Let  Q  be  the  number  of  samples  representing  h(t). 

2.  Choose  N  according  to  the  table  in  Figure  16. 

3.  Form  the  sampled  periodic  function  h(k) : 

h(k)  =  h(kT)  k  =  0,1,...,   Q  -  1 

=  0  k  =  Q,   Q+i N-  1 

4.  Compute  the  discrete  Fourier  transform  of  h(k) . 

N-i 

H(n)  =  E  h(k)  exp(-j2Trnk/N) 
k=0 

5.  Form  the  sampled  periodic  function: 

x±(k)  =  x(kT)     k  =  0,1,...,  N  -  Q 

=  0        k  =  N  -  Q  +  1,...,  N  -  1 

6.  Compute  the  discrete  Fourier  transform  of  x. (k) . 

N-l 
X.  (n)  =  E  x.  (k)  exp(-j2rrnk/N) 
1      k=0 

7.  Compute  the  product: 

Y±(n)  =  X1(n)H(n) 

8.  Compute  the  inverse  discrete  transform  of  Y. (n) , 

N-l  1 
yi(k)  =  E  jf  Y  (n)  #  exp(-j2nnk/N) 
n=0 

9.  Repeat  steps  5-8  until  all  sections  are  computed. 
10.  Combine  the  sectioned  results  by  the  relationships: 

y(k)  =  yi(k)  k  =  0,1, ...,N  -  Q 

y('fcW-QH)  =  y^k+N-0^-1)  +  y2(k)  k  =  0,l,...,N-Q 

y(k+2(N-Qfl))  =  y2(k+N-Q+l)  +  y3(k)  k=0,l,...,N-Q 
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